"""
Kuramoto-Sivashinsky PDE implementation

The Kuramoto-Sivashinsky equation describes chaotic spatiotemporal dynamics:
    ∂u/∂t = -ν∇⁴u - ∇²u - (1/2)|∇u|² + η

This is a fourth-order nonlinear PDE that can exhibit chaotic behavior.
"""

using ..Fields
using ..Grids
using Random

"""
    KuramotoSivashinskyPDE <: AbstractPDE

Kuramoto-Sivashinsky equation for pattern formation and chaos.

The mathematical definition is:
    ∂u/∂t = -ν∇⁴u - ∇²u - (1/2)|∇u|² + η(r,t)

where u is the field variable, ν controls the strength of the fourth-order term,
and η is optional Gaussian white noise.

# Parameters
- `nu`: Fourth-order diffusion coefficient (default: 1.0)
- `noise`: Noise variance (default: 0.0)
- `bc`: Boundary conditions for u (default: "auto_periodic_neumann")
- `bc_lap`: Boundary conditions for ∇²u (default: same as bc)

# Properties
- Fourth-order nonlinear PDE
- Can exhibit spatio-temporal chaos
- Dissipative system with finite-dimensional attractor
- Scale-invariant in certain parameter regimes
"""
struct KuramotoSivashinskyPDE <: AbstractPDE
    nu::Float64
    noise::Float64
    bc::String
    bc_lap::String
    rng::Random.AbstractRNG
    
    function KuramotoSivashinskyPDE(; nu=1.0, noise=0.0, 
                                  bc="auto_periodic_neumann", bc_lap="", 
                                  rng=Random.default_rng())
        bc_lap_final = isempty(bc_lap) ? bc : bc_lap
        new(nu, noise, bc, bc_lap_final, rng)
    end
end

"""
    evolution_rate(pde::KuramotoSivashinskyPDE, state::ScalarField, t::Float64)

Compute the evolution rate for the Kuramoto-Sivashinsky equation.

The equation is: ∂u/∂t = -ν∇⁴u - ∇²u - (1/2)|∇u|² + η

Implementation strategy:
1. Compute ∇²u (second-order diffusion term)
2. Compute ∇⁴u = ∇²(∇²u) (fourth-order dissipation term)
3. Compute |∇u|² (nonlinear term)
4. Combine all terms with appropriate signs
"""
function evolution_rate(pde::KuramotoSivashinskyPDE, state::ScalarField, t::Float64)
    u = data(state)
    grid_obj = grid(state)
    
    # Step 1: Compute ∇²u
    u_laplace = laplacian(state)
    
    # Step 2: Compute ∇⁴u = ∇²(∇²u)
    # Apply Laplacian again to get fourth-order term
    u_biharmonic = laplacian(u_laplace)
    fourth_order_term = -pde.nu * u_biharmonic
    
    # Step 3: Second-order diffusion term -∇²u
    second_order_term = -u_laplace
    
    # Step 4: Nonlinear gradient term -(1/2)|∇u|²
    grad_u = gradient(state)
    
    # Compute |∇u|² for each grid point
    grad_magnitude_sq = zeros(size(u))
    if isa(grad_u, VectorField)
        # For VectorField, sum squares of components
        grad_data = data(grad_u)
        if ndims(grad_data) >= 2  # Multi-dimensional gradients
            for i in 1:size(grad_data, 1)  # Loop over spatial dimensions
                grad_magnitude_sq .+= grad_data[i, :].^2
            end
        else
            grad_magnitude_sq = grad_data.^2
        end
    else
        # Handle case where gradient returns array of arrays or similar
        for component in grad_u
            grad_magnitude_sq .+= component.^2
        end
    end
    
    nonlinear_term = ScalarField(-grad_magnitude_sq / 2, grid_obj)
    
    # Step 5: Combine all deterministic terms
    result = fourth_order_term + second_order_term + nonlinear_term
    
    # Step 6: Add noise if specified
    if pde.noise > 0
        noise_data = sqrt(pde.noise) * randn(pde.rng, size(u))
        noise_term = ScalarField(noise_data, grid_obj)
        result = result + noise_term
    end
    
    return result
end

"""
    has_noise(pde::KuramotoSivashinskyPDE)

Check if the PDE includes stochastic noise terms.
"""
has_noise(pde::KuramotoSivashinskyPDE) = pde.noise > 0

"""
    is_chaotic(pde::KuramotoSivashinskyPDE)

The Kuramoto-Sivashinsky equation typically exhibits chaotic behavior.
"""
is_chaotic(pde::KuramotoSivashinskyPDE) = true

"""
    get_lyapunov_exponent_estimate(pde::KuramotoSivashinskyPDE, domain_size)

Rough estimate of the largest Lyapunov exponent based on domain size.
For KS equation: λ ≈ 0.08 - 0.1 depending on parameters and domain.
"""
function get_lyapunov_exponent_estimate(pde::KuramotoSivashinskyPDE, domain_size)
    # Empirical estimate - actual value depends on specific parameters
    base_exponent = 0.09
    # Larger domains tend to have slightly larger Lyapunov exponents
    size_factor = log(domain_size / (2π)) / 10
    return base_exponent + size_factor
end

"""
    get_energy(pde::KuramotoSivashinskyPDE, state::ScalarField)

Compute the kinetic energy: E = (1/2)∫u² dV
"""
function get_energy(pde::KuramotoSivashinskyPDE, state::ScalarField)
    u = data(state)
    return sum(u.^2) / (2 * length(u))
end

"""
    get_dissipation_rate(pde::KuramotoSivashinskyPDE, state::ScalarField)

Compute the energy dissipation rate based on the current state.
"""
function get_dissipation_rate(pde::KuramotoSivashinskyPDE, state::ScalarField)
    u = data(state)
    
    # ∇²u contribution
    u_laplace = laplacian(state)
    laplace_dissipation = sum(data(u_laplace).^2)
    
    # ∇⁴u contribution (bi-harmonic dissipation)
    u_biharmonic = laplacian(u_laplace)
    biharmonic_dissipation = pde.nu * sum(u .* data(u_biharmonic))
    
    return laplace_dissipation + biharmonic_dissipation
end

"""
    is_stiff(pde::KuramotoSivashinskyPDE)

The KS equation is stiff due to fourth-order derivatives.
"""
is_stiff(pde::KuramotoSivashinskyPDE) = true

"""
    requires_high_accuracy(pde::KuramotoSivashinskyPDE)

KS equation requires high accuracy to capture chaotic dynamics correctly.
"""
requires_high_accuracy(pde::KuramotoSivashinskyPDE) = true

"""
    get_time_scale(pde::KuramotoSivashinskyPDE, grid_spacing)

Estimate the characteristic time scale for the KS equation.
The fourth-order term dominates the time scale.
"""
function get_time_scale(pde::KuramotoSivashinskyPDE, grid_spacing)
    # Time scale is roughly L⁴/ν where L is grid spacing
    return grid_spacing^4 / pde.nu
end

"""
    expression(pde::KuramotoSivashinskyPDE)

Return a string representation of the KS PDE.
"""
function expression(pde::KuramotoSivashinskyPDE)
    terms = String[]
    
    # Fourth-order term
    if pde.nu != 0
        nu_str = pde.nu == 1.0 ? "" : "$(pde.nu)"
        push!(terms, "-$(nu_str)∇⁴u")
    end
    
    # Second-order term
    push!(terms, "-∇²u")
    
    # Nonlinear term
    push!(terms, "-0.5|∇u|²")
    
    # Noise term
    if pde.noise > 0
        push!(terms, "√$(pde.noise)η")
    end
    
    return join(terms, " ")
end

# Export the PDE type and analysis functions
export KuramotoSivashinskyPDE, has_noise, is_chaotic, get_lyapunov_exponent_estimate,
       get_energy, get_dissipation_rate, get_time_scale